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This paper is devoted to the development of a theoretical and computational framework to 
efficiently sample the statistically significant thermally activated reaction pathways, in multi- 
dimensional systems obeying Langevin dynamics. We show how to obtain the set of most prob- 
able reaction pathways and compute the corrections due to quadratic thermal fluctuations around 
such trajectories. We discuss how to obtain predictions for the evolution of arbitrary observables 
and how to generate conformations which are representative of the transition state ensemble. We 
present an illustrative implementation of our method by studying the diffusion of a point particle 
in a 2-dimensional funneled external potential. 



I. INTRODUCTION 

The problem of characterizing thermally activated transitions in systems with a large number of degrees of freedom 
plays a pivotal role in the study of many physical phenomena. From a theoretical perspective, the most natural 
strategy to describe the dynamics of a generic system is to integrate numerically the equations of motion, e.g. via 
Molecular Dynamics (MD) simulations. In principle, once a statistically significant sample of MD transitions has been 
obtained, it is possible to derive predictions for the time evolution of the observables, which can be matched against 
experimental measurements. 

Unfortunately, if the system contains a large number of degrees of freedom, or if the interesting transitions are rare 
events, characterized by free energy barriers of more than few units of fcs T, MD approaches become extremely costly 
and often impracticable. This is essentially because most of the computational time is invested in describing the 
motion of the system while it is exploring the portions of configuration space associated to the stable and meta-stable 
states. On the other hand, one is mostly interested in the information encoded in the reaction pathways joining the 
. different states. 

A prototype of problem in which the MD approach is unfeasible is the study of the protein folding reaction. In 
^ this case, the typical number of atoms to be simulated is of order 10*, the microscopic time-scale associated with the 
motion of the torsional angles of the polypeptide chain is of the order of the ns, i.e. many orders of magnitude smaller 
than the macroscopic time-scale associated with the inverse folding rate, which is in the range ms — s. State-of-the-art 
all-atom MD simulations cover time intervals in the range 10 — 10^ ns (see e.g. [l| and references therein) and can be 
useful only in investigating the reactions leading to local secondary structures formation, i.e. a-helices and /3-sheets. 

Given the limitations of the available MD simulations, several alternative approaches have been recently developed, 
to deal with the problem of characterizing the entire reaction in a high-dimensional space, using available computers. 
A drawback of these methods is that they rely either on ad-hoc assumptions about the dynamics obeyed by the system 
or on a specific choice of reaction coordinates . 
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In this work, we develop a theoretical framework which we call Dominant Reaction Pathways (DRP), which allows 
to rigorously identify the statistically significant transition pathways, without relying on any choice of reaction coor- 
dinates. Most importantly, the DRP approach is based on the low temperature expansion of the Langevin dynamics 
and, therefore, it does not involve any uncontrolled ad hoc approximation. 

A major advantage of the DRP approach is that it avoids investing computational time in simulating the local 
thermal motion in the metastable configurations. This is possible because the key equation can be formulated in a 
form which does not depend explicitly on the time variable. As a result, the computational difficulties associated with 
the existence of different time scales in rare thermally activated reactions, are naturally and rigorously bypassed. 

Some of the ideas underlying the DRP method were first introduced in [1, Q , in the context of the study of the 
protein folding reaction. In those works, we showed that it is possible to efficiently determine the most probable 
conformational transition pathways, both in the case of a coarse-grained and of a fully atomistic representation of 
a molecule. We now generalize this method to a much larger class of noise-driven phenomena in high-dimensional 
systems, and we present several theoretical developments. In particular, we discuss how it is possible to obtain 
predictions for the evolution of arbitrary observables, while systematically including the effect of thermal fluctuations 
around the most probable reaction paths. With this important development, the method does no longer focus on just 
the most probable transition trajectories in configuration space, but rather it allows to sample the entire ensemble of 
statistically relevant reaction pathways. 

In the present paper we give a self-contained presentation of the theoretical formalism and of the relevant algorithms 
and we shall illustrate how these can be implemented in practice, by studying the simple and intuitive problem of 
the diffusion of a point-particle in a 2-dimensional external potential. In a forthcoming publication, we shall present 
and test the application of the DRP method to the study of noise-driven conformational transitions of a polypeptide 
chain. 

This paper is organized as follows. In section |TT] we review the path-integral formulation of stochastic Langevin 
diffusion and in lllll we relate it to the problem of computing the time evolution of experimentally observable quantities. 
In section IIVI we provide a presentation of low-temperature expansion underlying the DRP approach and provide 
a rigorous definition of the ensemble of statistically significant reaction pathways. Sections |V] and IVII focus on the 
problem of efficiently sampling the ensemble of statistically significant reaction paths and predicting the time evolution 
of the observables. The subsequent section IVTIl is devoted to the problem of generating an ensemble of configurations 
which is representative of the transition state ensemble. In section lVlIII we provide a simple example of the application 
of the DRP method to the problem of the diffusion in a 2-dimensional funneled potential landscape. The main results 
and conclusions are summarized in llXI 



II. PATH INTEGRAL FORMULATION OF THE LANGEVIN DYNAMICS 



The starting point of our discussion of thermally activated transitions is the Langevin equation, which describes 
the out-of-equilibrium dynamics of a system coupled to a thermalized heat bath, 

mx = --/x-VU{x)+({t), (1) 

where x = (xi, X2, xn) is in general a vector specifying a point in configuration space (for example, the coordinates 
of all the atoms in a molecule), U{x) is a the potential energy and £,{t) is a memory-less Gaussian noise with zero 
average. For sake of simplicity and without loss of generality, in this section we present the formalism in the case of 
one particle in one spacial dimension. The generalization to the multi-dimensional case is straightforward. 

The term on the left hand side of Eq. ^ introduces incrtial effects which are damped on a time-scale m/^ = td. 
In the following, we shall focus on physical systems for which the damping time scale tjj is much smaller than the 
microscopic time-scale associated with the local conformational motion. For example, in the case of the dynamics of 
a typical amino-acid chain, is of order of the fraction of a ps, hence much smaller than the relevant microscopic 
time-scale associated to the dynamics of the torsional angles, which takes place at the ns time-scale. 

For t':^ td the Langevin Eq. ((T|) reduces to 

| = -^^Vf/(.)+,W, (2) 
where D = is the diffusion coefficient which we shall assume to be independent on the particles position^ and 

m-v ^ r- 



^ For the sake of simplicity, we shall further assume that all the particles in the system have the same diffusion coefficient. The general- 
ization to models in which each constituent has a different diffusion coefficient is trivial. 
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r]{t) = i ^(t) is a Gaussian noise that satisfy the fluctuation-dissipation relation, 

{7jitMt')) = 2D6it-t'), (3) 
The probability distribution sampled by the stochastic differential Eq. ([2]) satisfies then the Fokker-Planck equation 

4- P{x, t) = DV (-^ V U{x) P{x, t))+DV' P{x, t). (4) 



dt ' ' ' \kBT 

A universal property of all the solutions of the Eq. ([4]) is that, in the long time limit, they converge to the Boltzmann 
weight, regardless of their initial conditions, 

P(x,t) —> e ^ \ (^5j 

Hence, the Fokker-Planck equation accounts for the approach to thermal equilibrium with the heat-bath. By per- 
forming the substitution 

P{x, t) = 6"^^^"=^ t) (6) 
the Fokker-Planck Eq. Q can be recast in the form of a Schrodinger Equation in imaginary time: 

-^^(x,t) = 7?e//^(a;,i), (7) 

where the effective " Quantum Hamiltonian" operator reads 

Heff - -D\I^ + Veff{x), (8) 

and the effective potential is defined as 

= [{^U{x)f-2kBTV^U{x)\ (9) 

Hence, the problem of studying the diffusion of a classical particle with diffusion constant £) in a heat-bath at 
temperature T can be mapped into the problem of determining the quantum-mechanical propagation in imaginary 
time of a virtual system, defined by the effective quantmn Hamiltonian ([5]), interacting with the effective potential 

Based on the analogy with quantum mechanics, it is immediate to obtain a path-integral representation of the 
solution of subject to the boundary conditions x{—t/2) = Xi and x{t/2) = xj: 



Pixf,- 
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{U{x,)-U{x,}) r' ^ - f% dr (4if^ + Ve//[^(r)]) 



Vx{t) e y ^ "J. (10) 



This equation expresses the fact that the conditional Fokker-Planck probability is formally equivalent to a quantum- 
mechanical propagator in imaginary time. 

In order to illustrate the physical content of the effective potential Veff{x) defined in Eq.®, let us consider the 
probability for the system to remain in the same configuration x, after an infinitesimal time-interval dt. One finds 

P(x,di|x,0) - e-^^^^^^)"**. (11) 

Hence, the function Veff{x) measures the tendency of Langevin trajectories to evolve away from a given configuration 
X. In particular, the configurations where the effective potential is smallest are those where a system will be trapped 
for the longest time. Conversely, the residence time will be shortest at configurations with the largest Notice 
that such an effective interaction diverges in the low-temperature limit. Such a property of the effective potential 
allows to interpret physically the two terms in Eq. ([9]). The first term is always positive and provides the tendency of 
the system to move away from regions of configuration space where the forces are very large. In the low-temperature 
limit this is the driving effect. On the other hand, at higher temperatures, the interaction is dominated by the second 
term in ([5]), which reduces the value of Veff in the vicinity of narrow minima of the physical potential U, i.e. when 



FIG. 1: Left panel; an potential energywith two degenerate minima. Right panel: the corresponding effective potential 14// (a;). 



V^C/ is large and positive. This implies that the system will be trapped for long times in regions of configuration 
space characterized by a small conformational entropy. 

These observations about the physical interpretation of Veff{x) are illustrated by the example presented in Fig. [1] 
The left panel displays a potential function U {x) characterized by two degenerate minima, the one on the left hand 
side being broader than that on the right hand side. The right panel displays the corresponding effective potential 
(obtained by setting both ksT and D to 1). The configurations in the bottom of the broad minimum have a larger 
tendency to diffuse away, in an infinitesimal time, because there is a large number of nearby configurations which are 
thermally accessible. On the other hand, the configurations which are thermally accessible from the narrow minimum 
on the right are very few, because the potential raises very rapidly. As a result, the effective potential Veff{x) does 
not display degenerate minima, rather it has a pronounced global minimum in the proximity of the right minimum 
of U{x). Note also that, the effective potential displays maxima where the forces are largest. 



III. AVERAGE TIME EVOLUTION OF OBSERVABLES 



In many reactions of interest, experiments cannot probe the kinetics of a single transition, but only measure the 
average of a given quantity over a large number of independent transitions. For example, fluorescence experiments 
on the kinetics of protein folding measure the time evolution of the average relative distance between two specific 
amino-acid in the chain, during the reaction — see e.g. and references therein — . If the system relaxes to thermal 
equilibrium, such quantities converge to their average value in the canonical ensemble. In this section we discuss how 
to relate the average time evolution of an observable to the path integral representation of the conditional Fokker- 
Planck probability . We do so by further exploiting the formal connection between the stochastic and the quantum 
dynamics. 

Let us begin by introducing an Hilbert space Ti. spanned by the eigenvector of configuration operator 

X\x)=x\x). (12) 

Here x denotes a point in the configuration space. Let 0{x) be the real-valued observable of interest, which we will 
assume to be a function of the system configuration only. For example, 0{x) could be the position of a particle or 
the distance between two monomers in an amino-acid chain. We can associate to 0{x) an operator O, acting on the 
Hilbert space Ti: 

d\x) ^ 0{x)\x), (13) 
We are now in a condition to define the time-dependent average of the observable {0{t)): 

J dxe >'bt (xle-^^ff^xo) 

where t' is a arbitrary instant in the interval — 1/2 < t' < t/2. In practice, one is mostly interested in transitions 
which can take an arbitrary long time, i.e. those for which the total time interval diverges, t oo. 
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A more physically transparent representation of (jM]) is obtained by inserting a resolution of the identity, in terms 
of the eigenvalues of the configuration operator, \x'): 

Lf(x)-[/(xn)-[/(x') + L'(x') ' , , ,^ » - , / , , 

_ Jdxjdx'e~ (a:|e-^=//('/2-*) |2;')(a;'|Oe-^-"(*+*/2)|xo) 

\0[t )) — 



Jdx P{x,t/2\xQ,-t/2) 
J dx' J dxP{x,t/2\x' ,t')0{x')P{x' ,t'\xo, ~t/2) 
Jdx P{x,t/2\xo,-t/2) ■ 



(15) 



It is easy to check that this out-of-equilibrium average definition yields the correct result, in the small and large 
time limit. In fact, if the average is evaluated at initial time — i.e. for t' {—t/2) — it returns the value of the 
observable in the initial configuration, i.e. 

^^]im^^{0{t')) = 0(.To). (16) 

On the other hand, in appendix [Xl we prove that if the average (jl5l) is evaluated in the long time limit — i.e. for 
t' — > {t/2) and t —^ oo — it coincides with the equilibrium Boltzmann average, 

, lim (17) 
*'-*/2,*-oo Jdxe-^ 

The path integral representation of (jl5p is immediately derived from the path integral representation of 
P{x',t'\xo,-t/2) and P{x,t/2\x' ,t'), as in Eq. We find 



Idxe-^C^Vxjr) 0[x(t')] ^ 
~ Jdxe-^ Cvx{r) e-^-^>^ m^y^nMr)^) ' ^^'^ 

The discretized representation of (fT8|) is 



! dxe-^ Yf,-,' [J dx,] 0{x{)e 



m')) 



lim 



N 51]/ 



^dxe-^Ylt,' Udx,]e 



(19) 



where i is the discretized time step closest to t', i.e. t' ~ ij^. 

Monte Carlo simulations offer in principle a scheme to compute the average (|19p . If the interesting transitions are 
not rare events, then the number of uniform time discretization steps required to converge to the continuum limit is 
expected to be moderate. On the other hand, if the interesting transitions are rare, i.e. they occur at a time scale 
much longer than the microscopic time scale at which the local dynamics takes place, then the number of uniform 
discretization steps required to achieve the continuum limit would diverge exponentially. In this regime, the direct 
evaluation of the average ^8]) by Monte Carlo approaches rapidly becomes unfeasible. 

In the following sections we shall discuss how, in some appropriate regime, the sampling of the most important 
reaction trajectories and the calculation of the evolution of the observables can be efficiently performed by applying 
the instanton calculus, a technique which was originally developed in the context of the functional integral description 
of quantum tunneling. 



IV. STATISTICALLY SIGNIFICANT REACTION PATHWAYS 



Owing to the stochastic nature of the underlying dynamics, thermally activated reactions always take place through 
different trajectories in configuration space. However, not all such trajectories can be considered as truly inequivalent 
reaction pathways. An ensemble made by an infinite number of paths will exhibit similar features. For example, 
during two transitions, the system may cross the same sequence of saddle points or, in the specific case of macro- 
molecules, may undergo the same pattern of contact formations. Clearly, in these cases, it makes sense to view 
such reaction trajectories as different realizations of the same transition. On the contrary, trajectories which cross 
different saddles, or conformational transitions characterized by different sequences of contact formations represent 
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truly distinct reaction pathways. In general, the distinct pathways will not all be equally probable and will therefore 
contribute differently to the reaction rate. 

These observations illustrate two main problems which emerge when attempting a characterization of the significant 
reaction pathways in high dimensional systems. The first consists in providing a rigorous notion of reaction paths, i.e. 
one which does not rely on an a priori choice of reaction coordinates. The second problem regards the identification 
of the ensemble of most statistically relevant pathways, without having first to generate a huge sample of reaction 
trajectories by direct integration of the equation of motion. Once such a path has been identified, it is possible to 
make prediction for the average evolution of the observables and compare with experiments. 

In this section, we show how the path integral formulation of Langevin dynamics provides a natural way of defining, 
determining and characterizing the ensemble of most statistically significant reaction pathways, and to compute the 
time evolution of arbitrary observables, including the effect of thermal fluctuations, without relying on any prior choice 
of reaction coordinate. 

We start from Eq. (|10p . which expresses the probability to make a transition from an initial configuration Xi 
to a final configuration xj, in terms of a sum over all possible transition pathways connecting these two points in 
configuration space. Since each path is assigned a well-defined probability, e~^''ff^^\ the paths with the highest 
probability are those of minimum effective action Seff, i.e. those satisfying the Euler-Lagrange equations of motion 



2D 

with boundary conditions 



-^£ = VK//(x), (20) 



Xf. (21) 



These equations describe the classical (Newtonian) dynamics of a set of particles in the potential —Veff{x). In the 
following, we shall refer to the solutions of ([20]) obeying (PT|) as the dominant reaction trajectories. They correspond 
to instantons in the language of quantum tunneling theory. We stress that, in order to correctly account for entropic 
and energetic effects it is crucial to include the Laplacian term in the effective potential. As it was clearly shown in @], 
approaches which do not account for the Laplacian contribution do not correctly identify the saddle-point pathways. 

For a given choice of the boundary conditions, Xi and Xf, there are in general several dominant trajectories. Along 
with the exact solutions of the equation of motion associated to the effective action Seff, there are also gwasi-exact 
solutions. For example, if the relevant transitions involve overcoming an energetic barrier, any combination of exact 
independent solutions of the equations of motion which cross the barrier back and forth several times represents itself 
a solution, up to terms which are linear in the temperature and exponentially small in the height of the barrier 
Such multiple-crossing trajectories correspond to multi instanton-antiinstanton solutions in the quantum mechanical 
language. 

Obviously, since the exact and quasi-exact dominant trajectories form a functional set of zero measure, the proba- 
bility to follow exactly any of them during a transition is zero. A significant contribution to the transition probability 
comes from the set of all the paths which lie in the functional vicinity of a given dominant trajectory and differ 
only by small thermal fiuctuations. These paths form an ensemble of thermally equivalent realizations of the same 
transition. By contrast, the statistical weight of the huge set of pathways which are not in the functional vicinity of a 
dominant trajectory is exponentially suppressed by the factor e~^^ff . Hence, such paths do not provide a significant 
contribution to the path integral (jlOp and are very unlikely to occur. 

Let us take a closer look on the parameter controlling the regime in which the contribution of the saddle-points in 
the path integral dominates over that of the stochastic fluctuations. We begin by observing that the Fokker-Planck 
Eq. is invariant under a simultaneous rescaling of the diffusion coefficient and of the time variable, 

D^aD t^ -t. (22) 
a 

Hence, without loss of generality, we can restrict our attention to the Fokker-Planck Eq. in the form 

|- Pix,t)=.ViVU{x) P{x,t))+0V^ P{x,t). (23) 
ot 

where ^ is a parameter proportional to ksT. The path integral representation of the Fokker-Planck equation in the 
form ^ is 
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Notice that in this expression the coefRcient 20 plays the role of the combination fi/m in the quantum mechanical 
path integral. Hence, in general, the contribution of the paths in the vicinity of a saddle-point becomes important at 
low temperatures, when the parameter 9 is small. This regime corresponds to the semi-classical limit, in the quantum 
language 3l • 

In the low temperature limit, it is possible to express the transition probability (|10[) as a sum over a countable set 
of contributions i = 1,2, each of which accounts for the pathways which are close to a corresponding dominant 
trajectory, denoted as Xi(r). To obtain such a representation, let us take a completely general transition pathway 
x{t) and re- write it as the sum of the \-th dominant trajectory Xi(T) and a fluctuation j/(r) around it: 

x(t)=x,{t)+v{t). (25) 

Note that, by construction, the fluctuation function ?;(t) satisfies the boundary conditions y(T/2) — y(— t/2) = 0. 

We can now represent the contribution to the probability ()10p of the paths near the i—th dominant trajectory in 
the saddle-point approximation, i.e. by functionally expanding the action around a;i(r): 

1 /-t/Z pt/2 c2o r^l 



1 



2 

where we have introduced the operator F, defined as 

S^Seffix 



^e//[5.] + ^ /*^' dryir) F[x] ^(t), (26) 

t/2 



dxiT')SxiT") 



-D-^ + ^effin Kr - r') ^ F[n Sir r'), (27) 



which governs the dynamics of the fluctuations around Xi. Note that the expansion of the action does not contain 
linear terms in the fluctuation field ^(t), because the path around which we expand is a solution of the equations of 
motion. 

We now expand the fluctuation function y(r) in the basis of the (real) complete set of eigenfunctions of F: 

yir) = ^c„x„(t), F\xi]xn{T) = A„a;„(T). (28) 

n 

where the eigenfunctions Xn are chosen so as to satisfy the boundary conditions a;„(±i/2) = and the normalization 
condition 

t/2 

dTXra{T)Xn{T) = 5mn (29) 

-t/2 

To second order in the thermal fluctuations, the action reads 

Seff[x\= Seff\x,]^Yllcn\l + ... (30) 

n 

and the measure of the functional integral can be re-written as 

With this choice, the contribution of the thermal fluctuations around the i ^th dominant trajectory to the transition 
probability can be finally written as 

P.{xs,t/2\x.,t/2) ^ e-^(^(^^)-^(^'» e-^=/^[^'l H TT -J^e"^-^" 

■/-oo „ \/{2itY 

e-S.„[^.]. (32) 



detF\xi\ 
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The formal expression ((32|) represents a reliable estimate of the contribution of the i — th dominant trajectory to 
the transition probability pH]) only if the total time interval t is chosen to be of the same order magnitude as the 
microscopic time associated to the local motion of the relevant degrees of freedom x. 

On the other hand, in the study of thermally activated reactions, we are interested in transitions which take place on 
a time scale much larger than the microscopic time scales. Under such circumstances, the prefactor in the expression 
([32|) diverges and the approximate expression in Eq. ((32|) is no longer accurate, and has to be corrected [§]. 

To see how this problem emerges, it is useful to analyze the behavior of p2p in the limit in which the time interval 
for the transition diverges t oo. In such a limit, the time-translation invariance of the solutions of Eq. (j20p is 
restored, hence denoting by to the time at which a transition occurs one has 



Let us now consider the function 



Xi{t,to) = Xi{t~to). (33) 
MT-to)^ ^ .. M-T-to) (34) 

2DSeff[x,\ 

It is immediate to verify that this solution is an eigenstate of the fluctuation operator with zero eigenvalue (zero-mode) . 
To this end it is sufficient to take the derivative of the equations of motion: 

i!:^-2DVeff{x^)=Q (35) 

which leads to 

= {-l/2Dd^, + Veff)xoiT) = ~F[x^]xo (36) 

The physical origin of this zero- mode resides in the fact that, in the large t limit, the functional minimum of the 
action does not reduce to a single path; indeed, there is an infinite family of degenerate paths. Physically, such a 
large degeneracy reflects the fact that the crossing of the barrier can take place at any time between — oo and oo and 
that single-crossing events can occur with equal probability at any time. 

In the presence of a zero mode, the integral over cq in (j32p is divergent and cannot be performed in the Gaussian 
approximation. This problem has been circumvented by replacing the integral over cq with an integral over the time 
at which the transition takes place (see e.g. [HI). It is easy to show (see Appendix [5| that the integral over dco can 
be related to the time at which the transition takes place. 



dco^ ,^2DSeff[x,]dto. (37) 
Thus, the integration over dcQ can be performed exactly, leading to our final result 

P^{xf,t/2\x^,t/2)^t A, e"5E^(^("^)-^("'» e-^=^^[^'J, (38) 

where 



Y det'^^[a;,] ^ ^ 

and the symbol det' F denotes the product of all strictly positive eigenvalues of the fluctuation operator. 

The complete representation of the transition probability P{xf,t/2\xi,t/2) in the saddle-point approximation is 
obtained by summing up the contributions associated to the functional vicinity of each minimum of the Seff action, 
i.e. 

P{xf,t/2\x^,-t/2) ~ ^P,(a;/, V2|x„-t/2). (40) 



We stress the fact that each of the Pi in Eq. ([40]) includes the contribution of an ensemble of reaction paths, which 
are close to the same dominant reaction path. Hence, in the regime in which the saddle-point approximation is reliable, 
each dominant reaction path can be considered as the representative element of an infinite set of reaction pathways, 
which are equivalent up to small thermal fluctuations. Different dominant trajectories represent inequivalent reaction 
pathways and the set of all such inequivalent paths characterizes the ensemble of all possible distinct and statistically 
relevant pathways through which the transition can occur. To lowest-order in the saddle-point approximation, an 
estimate for the relative probability of each class of reaction pathway is provided by the weight cxp[— S'^y^]. 
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In the remaining part of this section, we discuss how to use the saddle-point approximation to calculate the average 
time-evolution of an arbitrary observable, over a statistically significant ensemble of independent transitions. 

According to the definition ([T5|l . one needs to average over the contribution of the dominant reaction pathways 

connecting the same initial configuration .tq to an ensemble of final configurations distributed according to e ^'--b^. 
To lowest-order, it is immediate to obtain the saddle-point expression for {0(t')), as it amounts to retaining only the 
contribution of the dominant pathways: 



J dxe '"'Bi- J2^ 



e 



where the sum runs over the dominant paths ^^(t), which connect a;o to some a; in a time interval t. 

At the next-to-leading order in the saddle-point approximation, each exponential weight in the sum Ij4ip is multiplied 
by the pre- factors Ai, representing the measure of the space of thermal fluctuations around each saddle-point, 

^ Jdxe-^ J:^Omt')] A, e-^e//[^-^-'-°-*l ^ ^^2) 
j dxe~^^ yl- g-Se//[3r,;x,xo,t] 

In practice, computing explicitly the Ai coefficients is a very challenging task, except for very simple systems (see e.g. 
0). In addition, if the initial configuration xq is not completely determined, one needs to perform an average over 
its probability distribution. 

The next sections are devoted to the problem of determining numerically the ensemble of most statistically significant 
transition pathways and calculating the average time-evolution of reaction coordinates, using numerical techniques. 

V. NUMERICAL DETERMINATION OF THE DOMINANT REACTION PATHWAYS 

In order to see how the dominant reaction pathways can be numerically determined, let us begin by recalling that 
each of them is a solution of the equation of motion 

-L^(r) = Vye//(x). (43) 
The solutions of (|43| conserve an effective "energy" Eeff defined as 

£^e// = ^s'-K//(S). (44) 

Note that £^e// has the dimension of a rate and therefore does not have the physical interpretation of a mechanical 
energy. 

As a consequence of the conservation of E^ff^ the symplectic action ^e// evaluated along a dominant trajectory 
can be rewritten as 

SeffiXf,Xi,t) = -E{t)t + SH.j{Xf,Xi\Eeff{t)), (45) 

where SHjixf,Xi]Eeff) is the Hamilton- Jacobi (HJ) action defined as 

SHj{.Xf,x,-E,ff{i)) = ^ dl v'i?e//(i) + K//[5(0], (46) 

where dl = Vdx^ is a measure of the distance covered by each constituent during the transition. This action can be 
seen to be identical to the argument of the exponent in the wave-function, in the semi-classical WKB approximation 
of quantum mechanics. 

We stress the fact that, in general, the numerical value of the effective energy E^fj depends on the choice of total 
time of the transition time through the relationship 

''f 1 
dl , (47) 

v^DjEijjw+vijjwm 
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Eq. (|46p expresses the fact that the dominant trajectory in configuration space leaving from Xi and reaching Xf is 
a minimum of the HJ functional. Note that, at finite temperature, the minimum of the physical potential U{x) does 
not in general coincide with the minimum of the effective potential Veff{x). However, the displacement between the 
minimum of U{x) and of Veff{x) are of order fcsT. 

In we have shown that, in the context of the study of configurational transitions of macro- molecules, the 
possibility of adopting the HJ formulation of the dynamics leads to an impressive computational simplification of the 
problem. In fact, the total distance between two different molecular configurations is at most 1-2 orders of magnitude 
larger than the most microscopic length scale, i.e. the typical monomer (or atom) size. As a consequence, ^ 50 — 100 
discretized displacement steps are usually sufficient for achieving convergence in the discretization of the trajectories. 
This number should be compared with 10^^ time-steps which would be required to characterize a protein folding 
trajectory adopting the time-dependent formulation. 

Once a dominant trajectory x{l) has been determined, one can recover the information on the time evolution of the 
system during the transition by means of the well-known relationship 

r(x) ^ f dl ^ (48) 



where t{x) is the instant between —t/2 and t/2 at which the path visits the configuration x. Hence, through Eq. 
it is possible to identify the configurations which are most probabily visited at each instant of time, i.e. the function 
x{t). 

Let us now discuss the implications of the choice of the effective energy, Eeff. From Eq. (|48p it follows immediately 
that large effective energies are associated to short transition times. In addition, from Eq. l(46|l we see that in the 
infinite effective energy limit, trajectories become insensitive on the structure of the interaction and the dominant 
path is the straight line in the multi-dimensional configuration space. 

In thermally activated processes, we are interested in transitions which take place on a very long time scale, which 
decouples from the microscopic time scale associated to the local motion in the stable and meta-stable states. In a 
protein folding transition, the final state x{t) should be taken as the native state of the protein xj, that is the absolute 
minimum of U{x), or possibly a minimum of Veff{x) close to it. In order for the protein to stay in this state as long 
as possible, the velocity of the system at that point should be x{t) — Xf ~ 0, which implies 

^^ff = - ^^//(^(^)) = -^«//(^(^))- (49) 

In conclusion, we have found that in thermally activated rare transitions, the most probable reaction pathways are 
those which minimize the functional 

SH.jixf,Xi;-Veff[Xf]) = J dl^Veff[x{l)] - Veffixf], (50) 

In appendix [C] we present a simple analytic example of the evaluation of the estimate for the transition probability, 
to lowest order in the saddle-point approximation, using the HJ formalism. 

The numerical determination of the dominant reaction pathways involves searching for the minima of the discretized 
version of HJ functional (|46p . For example, if we are interested in molecular transitions, the discretized HJ functional 

N-l 

E V (^e//(n) - K//(l))AZ„,„+i + AP, (51) 

n 

where N is the total number of discretization steps and 



reads 



D ^ 

^ ^ i—l 



N \ 

EVj{/(xi(n),...,XAr(n)) - 2kBTj2^PiMn): -^^Nin)) 



(52) 

P = E(^k"+i-(^0)^ (53) 

n 

(A/)^„+i = Y,iMn+l)-Mn)r- (54) 
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Np is the number of particles (atoms or residues) in the molecule, A ln,n+i is the Euclidean measure of the n — th 
elementary path-step and A P is a penalty function which keeps all the length elements close to their average [2[. 
Although this term becomes irrelevant in the continuum limit, it plays an important role if the transitions are rare 
events. In fact, it enforces that the sampling is made at constant displacement steps rather than at constant time 
steps. The minimization of the discretized HJ functional (j51[) can be performed for example by means of a simulated 
annealing algorithm. 

We conclude this section by noting that the coordinate I appearing in can be interpreted as a reaction- 
coordinate. However, we stress that this quantity is not chosen a priori, but rather emerges naturally and in a 
self-consistent way from the path-integral formulation of the Langevin dynamics. The dominant path encodes some 
information about the structure of the free energy landscape along such a coordinate. In fact, one expects that the 
residence time in a configuration visited by a dominant path to be largest (smallest) in regions where the free-energy is 
smallest (largest). Indeed, numerical simulations of the conformational transitions of alanine dipeptide fully supported 
this expectation and showed that this way it was possible to locate very accurately the maxima (minima) of the 
free energy landscape. 



VI. THERMAL FLUCTUATIONS AROUND THE DOMINANT REACTION PATHWAYS 



In the previous section we have shown how to determine the set of dominant reaction pathways corresponding to 
a given choice of boundary conditions. In this section, we show how it is possible to include the effect of quadratic 
thermal fluctuations around a dominant path, by means of the Monte Carlo algorithm. For the sake of definiteness, we 
focus on molecular transitions. The method, however, holds in general for thermally activated reactions in arbitrary 
systems obeying Langevin dynamics. 

Let us denote by 

x{n) = (xi(n), ...,X7Vp(n)), (55) 

with n = 1, ..,iV, the dominant trajectory corresponding to a sequence of configurations of the molecule, and deter- 
mined by minimizing numerically the HJ action (|5ip . 

The time at which each configuration x{n) is visited during the transition can be obtained by computing the set of 
time intervals separating each of the path steps x{n) from x(n + 1): 

Ar„ „ I 1 = — = (56) 

' + ^AD{V,jf{x{n))-V,ff{x{l))) 

For example, the 10 — th configuration a;(10) is visited at the time r(10) = X^n*^^ Ar„,„+i. 

The sequence of time intervals obtained this way can be used to write a discretized version of the transition 
probability pO)) : 



P{xf,x^;T{N)) ^ / _[_[ [dxi(n)...dxA,p(n)J e L ^ "•"+^ ^ J (57) 

Eq. ([57)) represents a particularly convenient Trotter discretization of the original path integral pI7|) . In fact, the 
time intervals are chosen large, when the most probable trajectory at that time is evolving slowly. On the other 
hand, the time intervals become very small when the configurations in the dominant path are moving fast. Such 
a construction is reliable only in the low temperature limit, where the information encoded in the saddle-point is 
relevant. In this regime, the number of Trotter discretization steps N is the same as the number of path discretization 
steps in (|5ip and is much smaller than the number of equal time-discretization steps, which would be required if 
one had no prior knowledge of the relevant transition. From the computational point of view, the sampling of the 
ensemble of fluctuations around each saddle-point in ([57]) with a Monte Carlo algorithm is essentially as expensive as 
finding the dominant path itself. In fact, in both case one has to sample a space with iNp x N degrees of freedom. 

In order to single out the next-to-leading corrections in the saddle-point approximation of the propagator, one 
should in principle expand the action to quadratic order around the saddle-point first and then perform the MC 
sampling. However, in the regime of validity of the saddle-point approximation, the contribution to the propagator 
of the cubic, quartic and higher terms in the expansion of the action are small and lead to corrections which come in 
at next-to-next-to-leading order. Hence, in such a regime, one does not need to perform the expansion of the action. 
On the other hand, if this condition is not satisfied, i.e. if thermal fluctuations are large, then the entire saddle-point 
approach breaks down. The approximation scheme presented in this section can be applied also to derive estimate 
for the average value of observable, at next-to-leading order. 
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FIG. 2: Two dimensional energy landscapes considered in section IVIIII 



VII. CHARACTERIZATION OF THE TRANSITION STATE ENSEMBLE 

The DRP formalism provides a prescription to identify a set of conformations which can be taken as the most 
representative of the Transition State Ensemble (TSE). 

To see how is it possible to extract such information from the dominant reaction pathways, we adopt the definition 
TSE in terms of commitment analysis [isj. i.e. as the ensemble of conformations from which the probability to 
evolve into the native state equals the probability of collapsing into a denatured conformation. Once a dominant 
path has been determined, it is immediate to identify the conformation xts such that the probability in the saddle- 
point approximation to diffuse back to the initial configuration Xi equates the probability of evolving toward the final 
configuration Xf, i.e. 

P{xN\xts)^P{xD\xts)- (58) 

To lowest order in the saddle-point approximation, this condition leads to the simple equation: 

U{x,) - U{xf) 



2kBT 



= SHj{[x\\xts-,Xf) - SH.j{[x\;xts-,Xi). (59) 



The determination of the representative elements of the TSE can be generalized beyond the lowest-order, of the 
saddle-point approximation. To this end, we recall that the dominant reaction pathway is obtained in practice by 
minimizing numerically the discretized version of the HJ action (jSip . Hence, the configuration xts which is most 
representative of the TSE corresponds to a given path step, rits- Once the most significant thermal fluctuations 
around the dominant pathway have been determined, using the prescription described in section IVTl it is possible to 
determine a full ensemble of representative configurations of the transition state, 

TSEd {xi{nts), ■■■,XN{nts)}, (60) 

where the index l..iV labels all the fluctuation paths, and Xi{n) with i = 1, ..N are the paths determined with Monte 
Carlo sampling of the integral (|57|) . 



VIII. AN ILLUSTRATIVE EXAMPLE: NOISE-DRIVEN ESCAPE FROM A 2-DIMENSIONAL FUNNEL 

In view of applications of the DRP method to the study of transitions which take place in complex high-dimensional 
spaces, it is instructive to first illustrate how this approach can be implemented, in a very simple and intuitive case. 
In such a perspective, in this section we consider a 2-dimensional model, which displays some of the key features of 
the protein folding reaction. We consider a funneled energy function ?7(x), whose minimum is narrow and several 
ksT deep. Such a minimum is separated by a fiat circular region by an energetic barrier of some fcsT. The flat region 
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FIG. 3: A typical starting trajectory before and after applying the smearing algorithm. 



is corrugated by adding a hump or a hole, of few ksT. The system is confined by a symmetric harmonic potential, 
which becomes much larger that fc^T only at large distances from the center of the funnel. 
The specific functional form of the chosen potential energy function is 

C/(x) = C/^(x) + C/^(x) + i7^(x). (61) 
t/^(x) is the term providing the funnel with an energetic barrier at its border and reads 

[/^(x) = ae-V(2«^) _ 6e-'''/(2.= )^ (g2) 
where a — 45 kJmol^^, u — 25nm, b — 45 kJmol^^ and v = lOnm. U^{x) is the term providing the hump and reads 

UR{^)=ce-^^-^°^"/^^" (63) 

where Xg is the position of the corrugation with size w — 3nm and hight c = 12.5k J. 
(x) is a confining and symmetric harmonic trap, 

C/^(x) = ia>2(x2), (64) 

where lu = 0.11k Jmol^^nm'^ . The resulting energy landscape is presented in Fig. ^ 

With the protein folding problem in mind, we shall refer to the ensemble of configurations near the bottom of the 
funnel as to the native state. Conversely, the metastable fiat region between the harmonic trap and the energetic 
barrier near the entrance of the funnel will be referred to as the denatured state. Note that the conformational entropy 
of the native state is much smaller than that of the denatured state, as it is the case for proteins. 

We are interested in studying the thermally activated transitions connecting an initial denatured configuration 
Xi to a native conformation Xf. Due to presence of the hump, the landscape is not spherically symmetric and the 
structure of the relevant transition pathways depends on the temperature of the heat-bath. At asymptotically high 
temperatures, the diffusion becomes quasi-Brownian, as the reaction pathways are expected to be little sensitive to 
the local structure of the energy landscape. In such a temperature regime, the DRP approach is not expected to be 
reliable. Conversely, in the low temperature regime, the diffusion will be driven by the structure of the force field: the 
most statistically relevant transitions will avoid visiting the positive corrugation and prefer to spend longer fractions 
of the transition time in the bottom of the funnel. In this regime the bundles of paths around different dominant 
reaction pathways will have negligible overlap and the DRP method will be reliable. 

The main point of the DRP method is to first generate the set of most probable reaction trajectories and then to 
use such information to efficiently sample the thermal fluctuations around each of such trajectories. This way one 
generates an ensemble of transition pathways, whose high statistical weight is guaranteed a priori, by construction. 

Any algorithm for the minimization of the HJ functional has to be provided with a starting path connecting the 
initial and final points, from which it begins the search for the most probable trajectories. In order to sample as 
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efficiently as possible the space of transition pathways, one needs to generate a full ensemble of independent starting 
paths, each of which is defined in a relatively small number N of discretization steps. 

In this simple example, we generate such an ensemble of initial paths as follows. We run a set of MD simulations 
starting from the same denatured conformation Xj . The simulations are stopped when the system reaches the native 
state, i.e. when the potential energy takes a value which is less than 2kBT larger than than the potential at the 
minimum of the funnel^. 

The escape from the funnel represents thermally activated events, in which the microscopic time scale associated 
with the local motion in the native and denatured state is decoupled from the time scale associated to the inverse 
transition rate. Consequently, the initial trajectories obtained via MD simulations contain a very large number of 
discretization steps Nmd- On the other hand, the main advantage of the DRP method is that it allows to keep only 
a small number of path steps Njjfip, typically of the order of 50 — 100, for realistic molecular systems. Consequently, 
before beginning the minimization procedure, one needs to reduce the number of path steps in the initial MD trajectory 
from Nmd down to Njyfip ^ Nmd- 

In order to do so, we applied the following smearing algorithm to the initial path {x(i)}i=i...jvA,£, . At each iteration, 
we defined a new path {x'(i)}i=i.. jvMn : consisting of Nmd/'2 steps, in which each configuration was defined as average 
of two configurations of the input path at consecutive time steps, i.e. 

x'(^) = i(x(2i + l)+x(2i)). (65) 

After k of such decimation steps, one obtains a smeared trajectory consisting of Ndrp = 2^^Nmd steps. The shape 
of a typical initial trajectory before and after smearing is reported in Fig. [3l The smearing algorithm was used to 
reduce the number of steps in the trajectory from 0(10^) to 30 — 50. 

The minimization of the discretized HJ effective action was performed applying the following simple simulated 
annealing algorithm. After a preliminary thermalization phase, based on the Metropolis algorithm, we performed 10 
cooling cycles, which were stopped when the acceptance rate dropped below 10%. After each cycle we performed 
a re-heating cycle consisting of few Metropolis steps. At each step we updated the path globally, i.e. moving the 
position of the particle in each configuration of the path. At the end of each cooling cycle, the boldness of the Monte 
Carlo moves was adjusted, in order to keep the acceptance rate between 10% and 90%. The evolution of the HJ action 
during the annealing of one of the initial trajectories is reported in Fig. 2) For such a simple problem, convergence 
was typically achieved in few minutes on a single desktop. 

Once a dominant trajectory x{n), with n = 1, .., Nupp has been determined by minimizing numerically the HJ 
action (|5ip it is possible to sample the space of small thermal fluctuations around it by means of the Monte Carlo 
algorithm described in section IVIl i.e. by sampling the path integral ([57)) with a simple Monte Carlo simulation, 
using the dominant reaction pathway determined before as the starting configuration of the Markov chain. Given the 
simplicity of the system under consideration, it was sufficient to adopt a Metropolis algorithm. 

Let us now discuss the results concerning the diffusion from the denatured configuration Xi = (0, 50) [nrn] to the 
native conformation Xf — (0,0) [nm], in two different energy landscapes. In the first landscape, which we shall refer 




200 400 

Annealing steps 



FIG. 4: Typical evolution of the HJ action during the annealing algorithm. 



^ Such a prescription is suitable for the present study, aiming to illustrate how the DRP works. However, implementing the same procedure 
becomes very computationally demanding, when applied to generate starting paths in complicated multi-dimensional systems. In the 
discussion of the application of the DRP method to molecular systems, we shall adopt a different algorithm to generate the first path 
and choose a representative set of initial conformations Il4ll . 
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FIG. 5: Left Panel: dominant reaction pathway for the diffusion in the landscape A, obtained from 10 different starting 
trajectories at 300K (solid line) and 50A" (dashed lines). Right Panel: dominant reaction pathways for the diffusion in the 
landscape B obtained from 10 different trajectories at 50K (solid line) and 50K (dashed line). 



to as A, the hump in the flat region was removed and the diffusion took place in a smooth symmetric funnel. In the 
energy landscape B, a hump of bksT at T = SOOK was placed between the native and denatured conformations. The 
dominant reaction pathway for such a reaction have been found at different temperatures, T — 50K and T = 300K, 
by minimizing the HJ action using 10 independent initial trajectories. 

In the landscape A, the minimization algorithm always converged to the same dominant path, regardless of the 
temperature of the heat-bath and on the trajectory used as starting point of the annealing Markov chain — see the 
left panel of Fig. [5] — . Hence, for this transition, there is only one type of dominant path'^. 

In landscape B, the minimizations always converged to two different dominant reaction pathways: one passing on 
the right and one passing on the left of the barrier as shown in the right panel of Fig. [Sj Notice that the two dominant 
paths tend to converge, as the temperature is increased. This is well understood since, in the asymptotically large 
temperature limit, the structure of the landscape becomes irrelevant and the diffusion becomes purely Brownian. 

The temperature of the heat bath does not only affect the structure of the dominant reaction pathways. Most 
importantly, it governs the size of thermal fluctuations around the dominant pathways. The panels of Fig. [5] show the 
ensemble of statistically relevant fluctuations around the two dominant paths, obtained from the Monte Carlo sampling 




FIG. 6: Left Panel: Ensemble of statistical ffuctuations around the dominant reaction pathways at T = 300K. Right Panel: 
Ensemble of statistical fluctuations around the dominant reaction pathways at T = 50K. 



In the present analysis, we are not concerned with multi-instanton paths which escape from and return to the bottom of the funnel, 
several times. 
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FIG. 7: Left panel: energy landscape C. Right panel: residence time along the corresponding dominant reaction pathway. 

of the integral ([57)) . at temperature T = SOOif and T — 50i<r, respectively. As expected, at low temperatures, the 
fluctuations are small and the bundles of paths associated with the two dominant paths do not overlap. Consistently, 
the Monte Carlo sampling performed starting from one of the two dominant paths takes an exponentially large time 
to generate paths which cross the hump and explore the configurations in the other side of the valley. In this regime, 
it makes sense to distinguish between the two possible ways in which the transition can take place and the correct 
average trajectory is recovered only after one performs the average over the paths sampled around each dominant 
reaction pathways. 

On the other hand, as the temperature is increased and becomes of the order of the height of the barrier, the two 
sets of fluctuation trajectories begin to overlap and eventually become indistinguishable. This is evident from the 
second panel of Fig. [6l Notice that in all cases the DRP approach leads to the correct average trajectory which travels 
through the barrier in the middle. 

Let us now discuss the residence time in each of the configurations of a dominant reaction pathway. To this end, 
we consider a third landscape C which displays a negative depression, before the entrance of the funnel (left panel 
of Fig. [7]), corresponding to setting c = lOfcJ and w = Snm in Eq. ((63|) . In the right panel of Fig. [7l we plot 
the corresponding residence time along the dominant trajectory, as a function of the curvilinear abscissa which 
parametrizes the path. We clearly observe the presence of a local maxima of the residence time, when the particle 
travels near the depression. As expected, the residence time diverges in the native state, because of our choice of 
effective energy Eeff — ~Veff{xf). 

IX. CONCLUSIONS AND OUTLOOK 

In this work, we have provided a presentation of a theoretical framework to characterize the important thermally 
activated reaction pathways in multi-dimensional systems obeying Langevin dynamics. We have shown that, by 
applying the instanton calculus to the path integral formulation of the stochastic dynamics, it is possible to define, 
determine and characterize not only the most probable pathways, but also the ensemble of most statistically significant 
thermal fluctuations around them. 

From the theoretical point of view, our method has the advantage that it does not require any prior choice of 
reaction coordinate. However, a natural reaction coordinate, / emerges naturally from the HJ formulation of the 
dynamics of the dominant reaction pathways. In addition, this framework yields a simple procedure to generate sets 
of conformation which are representative of the transition state. 

From the computational point of view, the method has the advantage that it avoids investing CPU time in simulating 
the dynamics of the system in the metastable states, i.e. when it is not performing a transition to the final state. 

As a result, rather than simulating long Molecular Dynamics (MD) trajectories until a significant number of 
transitions occur, we suggest to first determine the most probable transitions and then to sample the most significant 
fluctuations around them. 
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APPENDIX A: LONG TIME LIMIT OF THE OUT-OF-EQUILIBRIUM AVERAGE 

Here we show that, if the average pS]) is evaluated in the long time limit — i.e. for t' {t/2) and t oo — it 
coincides with the equilibrium Boltzmann average, 

, li- (Al) 

We start from 

iJ(^)-iJ(xo) ^ . 

(0(0) = „:::„j-i"- 1-°' . (A2) 

We now insert a resolution of the identity in terms of a complete set of eigenstates of the effective quantum Hamiltonian 
In): 



U(a;)-l7(xo) 



S dx Y,n^~^"^ ^ ^'''''^ 1pnix)'lp„{xo) 

where we have introduced the real "wave functions" ipn{x) = {n\x) — ipn{x) = {x\n). 

In the long time limit, the contributions of all the eigenstates \n), except the ground-state are exponentially 
suppressed and one gets 

m)) = (A4) 

Jdxe ^"B-r" e--^o*V'(S(a;)f/'o(a;o) 
We now recall that for any systems obeying the Fokker-Planck Eq.Q, the ground-state of the effective Hamiltonian 

_£(xo) 

is V'o(a^) = e ^''^^ and the corresponding energy is null. This fact immediately yields Eq. p7j) . 



APPENDIX B: INTEGRATION IN THE ZERO-MODE SPACE 

In the presence of barrier crossing pathways, the WKB expression ([5^ breaks down and the integration over the 
zero mode coefficient cq cannot be performed in gaussian approximation, as 



ze 



7(2^ J-oo v/(2^ 



(Bl) 



Let us compare the integration the deviation dx from a given path which are caused by a small time translation 
of the time at which the transition takes place to ^ to + dto with the deviation caused by a small coefficient shift, 
Co —I- Co + dco- The shift in to leads to 



dxi 

dx — dto — —Xidto, (B2) 

dto 



while the shift cq leads to 



dx = ^dco = xodco = A / —-^ — x^dco (B3) 

dco y ^DSeff 

Equating the two variations and redefining the sign of co so that the two integrations have the same boundaries, we 
find 



dco = ^2DSeffdto. (B4) 



18 



APPENDIX C: A SIMPLE ANALYTIC EXAMPLE 



Here we illustrate the construction of the leading-order estimate of the Fokker-Planck conditional probability 

P{xft/2\x,i,—t/2) using the HJ formalism, in the simplest possible case, i.e. for pure Brownian motion. 

By rotational invariance, the Fokker-Planck probability is a function only of the distance L = \xf — Xi\ and of the 
time interval t. To lowest-order in the saddle-point approximation we have 

P{L,t) ~ e^''"'^*^*~^"-''^^''"'^'^^'^\ (CI) 

The dominant reaction pathway for this transition is obviously the straight line path, since the HJ action is proportional 
to the distance covered. Let us first evaluate the effective energy parameter as a function of the time interval , Eeff{t), 

// 1 r 

dl , = , (C2) 



which gives 



The HJ action is 



Eeffit) = (C3) 



SHj{E,ffit),L) = j^^ dlsj^{E,ff{t))=Lsj^^ 



We are therefore able to reconstruct the Fokker-Planck probability 

P{L,t) ~ const. ■ ero*~^ = const. ■ e"™* (C5) 

For such a simple problem the leading-order approximation yields the exact result. The normalization of the proba- 
bility is of course given by the quadratic fluctuations. 
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